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ABSTRACT 

A user’s guide for the program gmc.f is presented. The program is based on the generalized 
method of cells model (GMC), which is capable via a micromechanical analysis, of predicting 
the overall, inelastic behavior of unidirectional, multi-phases composites from the knowledge of 
the properties of the viscoplastic constituents. In particular, the program is sufficiently general to 
predict the response of undirectional composites having variable fiber shapes and different fiber 
arrays. 


INTRODUCTION 

The method of cells is a micromechanical model which has been shown to accurately 
predict the overall behavior of various types of composites from the knowledge of the consti- 
tuent properties. In particular, the method yields explicit effective constitutive equations for the 
inelestic behavior of metal matrix composites. The various capabilities of the method were 
demonstrated in a recent review paper and a monograph [1,2]. 

The overall behavior of inelastic, multi-phase, unidirectional fibrous composites generated 
by the method of cells from the knowledge of the properties of the individual constituents is 
displayed in terms of: 

• effective elastic moduli 

• effective coefficients of thermal expansion 

• effective thermal conductivities 

• effective stress-strain response in the inelastic region 
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Original method of cells 

In the original formulation of the method of cells, a continously-reinforced, unidirectional 
fibrous composite is modeled as a doubly-periodic array of fibers embedded in a matrix phase. 
Figure la. The periodic character of the assemblage allows one to identify a repeating unit cell 
that can be used as a building block to construct the entire composite. Figure lb. The properties 
of this repeating unit cell are thus representative of the properties of the entire assemblage. The 
unit cell consists of a single fiber subcell surrounded by three matrix subcells. Hence the name 
method of cells. The rectangular geometry of the repeating unit cell allows one to obtain an 
approximate solution for the stresses and strains in the individual subcells given some macros- 
copically homogeneous state of strain or stress applied to the composite. The approximate solu- 
tion to the thus posed boundary- value problem is, in turn, used to determine macroscopic (aver- 
age) or effective properties of the composite. 

The displacement field in the individual subcells is approximated in terms of a linear 

— (P) — (Y) 

expansion in the local coordinates x , x centered at the mid-point of a given subcell. 


u pY) =0) (^) + xf<t)^ ) + X3 Y VP Y) 


i = 1,2,3 


( 1 ) 


where cop are the displacement components at the center of each subcell (Py). Using the local 
strain-displacement equations for each subcell, one can easily show that the microvariables <f>P 
and \|/p are related to the subcell strains ep. For example. 



atop 

dxj 


eP = <|>P, ep = yp. 


( 2 ) 


with the shear strains given by similar expressions. 

The displacement expansion given by Equation (1) contains 9 unknown variables in each 
subcell of the repeating unit cell, i.e., three top ’s, three <J)p ’s and three \j/p Y) ’s, for a total of 
36 unknown variables per unit cell. These unknown variables are expressed in terms of the 
macroscopic or average strains applied to the composite by imposing the continuity of displace- 
ments between the individual subcells, as well as between the given repeating cell and the sur- 
rounding cells, together with the continuity of interfacial tractions. These continuity conditions 
are applied in an average sense due to the form of the diplacement representation in the indivi- 
dual subcells. Application of these continuity conditions leads to the following set of equations: 
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Continuity of interfacial displacements 


h] <j>P Y) + h 2 <t>P Y) = (hi + h2)-^- 


llV pi) + i 2V (p»= ( i 1+ i,) 


3o)i 

3xi 


- -, ( 


1 5(0j 8o>j 


2 dxj dx; 


) 


where ooj 11 ^ = 0)[ 12 ^ = cop 1 ^ - top*^ - C0j. 


Continuity of interfacial tractions 


_UY) -(2y) 

&2i - a 2i 

_(PD —(1*2) 

<*3i - <^3i 


( 3 ) 


(4) 

(5) 


Equations (3)-(5), together with the appropriate constitutive equations for the given subcell 
(Py), are used to express the microvariables in terms of the macroscopic strains, £jj, applied to 
the composite. The form of the constitutive equations used to describe the material behavior in 
the subcell (Py) is, 


o (pY) = C (pY) [e pY) - - a (PY) AT ] (6) 

(Py) __ (py) _p<Py) 

where a is the average stress, e is the average total strain, e is the average plastic 
strain, C ( ^ is the elastic stiffness matrix, a ( ^ Y - are the coefficients of thermal expansion, and AT 
is the temperature deviation of the material occupying subcell (Py). The volume averaged subcell 
stresses a^ Y) are defined in the usual way. 



+hp/2 +U/2 

-2- j } aW<K ? dxf 

V Py -hp/2 ~ly/2 


(7) 
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(where Vp 7 = hply), with similar definitions for the remaining averaged quantities. The volume 
averaged subcell stresses are obtained in terms of the microvariables in the following form: 



The above formulation is sufficiently general to admit any inelastic constitutive model for the 
average plastic strains in the subcell ((3y). Presently, the Bodner-Partom unified viscoplasti- 
city theory is used to model the inelastic effects. This constitutive model will be briefly 
described in the following section. 

The effective or average stress-strain equations for the composite are subsequently con- 
structed by applying the definition of the average composite stresses: 


where V = (h t + h2)0i + 12 ). The resulting effective stress-strain equations for the composite 
obtained from the above analysis have the form: 


where the elements of the effective elastic stiffness tensor C, the elements of the effective ther- 
mal expansion coefficient a. and the effective plastic strains e of the composite are expressed 
in terms of the elastic thermal properties and plastic strains of the constituent phases in closed 



( 8 ) 


a = C(e - e* 5 - aAT) 


( 9 ) 
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form. These expressions are lengthy and thus will not be reproduced here - they can be found in 
References [1,2]. 

Although the method of cells has been demonstrated in numerous experimental/analytical 
correlation studies to be an accurate and, at the same time, efficient tool for analyzing the inelas- 
tic response of metal matrix composites in a wide temperature range, it suffers from several 
drawbacks. The rectangular geometry limits the analysis of MMC’s to a limited number of fiber 
arrays. Only regular arrays with, at most, two different fiber spacings can be analyzed by the ori- 
ginal method. The use of four subcells precludes the possibility of inclusion of an interfacial 
region between the fiber and the matrix phases, as well as consideration of more complicated 
fiber shapes. The transversely isotropic behavior of unidirectional composites in the plane per- 
pendicular to the fiber direction is imposed artificially by setting hq = lj and Iq = b (see Figure 
1), and subsequently averaging out the effect of the square geometry of the unit cell. These limi- 
tations motivated the development of the generalized version of the method of cells described 
subsequently. 

Generalized method of cells 

In the generalized formulation, the repeating unit cell is subdivided into an arbitrary 
number of subcells. This generalization extends the modelling capability of the method of cells 
to include the following: 

• thermomechanical response of multi-phase, metal matrix composites 

• modelling of variable fiber shapes 

• analysis of different fiber arrays 

• modelling of porosities and damage 

• modelling of interfacial regions around inclusions, including interfacial degradation 

The generalization of the original method of cells allows approximate micromechanical 
analysis of more complicated periodic arrays such as the one shown in Figure 2. The repeating 
volume element used in the GMC micromechanical analysis is shown in Figure 3. It consists of 
NpxNy subcells. Each of these subcells is assumed to be occupied, in general, by an elastic- 
viscoplastic material. Thus the repeating volume element consists of NpxNy different viscoplas- 
tic materials, i.e., it represents a multi-phase, inelastic composite. 
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A given elastic-viscoplastic phase is characterized in the elastic region by five elastic con- 
stants: 


E a = axial Young’s modulus (sometimes called longitudinal or major) 
Ex = transverse Young’s modulus 
v A = axial Poisson’s ratio 
Vj = transverse Poisson’s ratio 
G a = axial shear modulus 


It is also characterized by two coefficients of thermal expansion: 


a A = axial CTE 
otj = transverse CTE 

The form of the constitutive equation used to described the material behavior in the subcell 
(( 3 y) in the generalized method of cells is precisely the same as for the original method of cells 
given by Equation (6). For a transversely isotropic material occupying the subcell ((fy), the form 
of the elastic stiffness matrix in Equation (6) is: 


C (PY) = 


Cn C12 C\2 0 0 0 

C12 C22 C23 0 0 0 


C12 C23 C22 0 

0 0 0 C44 

0 0 0 

0 0 


(Py) 


0 0 
0 0 
0 C55 0 


000c 


66 


( 10 ) 


where the X] direction is perpendicular to the plane of isotropy. The elastic stiffness matrix ele- 
ments are related to the more commonly used engineering constants in the following manner: 


C11 

= E A + 4kv 

C12 

= 2kv a 

C22 

= K 4- G T 

C23 

= K ~ Gj 

C44 

= G j 

C55 

= C 66 = G a 


where 
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and 


Gj = 


1 E t 

2 1 + Vj 


K = 


1 

4 l Ea 2 i 

[—(1 - v t)t v a ] 

2 Ex 


The form of the thermal expansion coefficient matrix is: 

a ( ^> = diag[aJA a? Y) , af# 0, 0, 0] (1 1) 

In the inelastic region the viscoplastic phase is represented at this time by the Bodner- 
Partom unified viscoplasticity theory. The local viscoplastic strain rate for a given inelastic sub- 
cell is expressed as 


• P (Py) a £ (Pt) 
e = A(py) s , p + Y * 2 


( 12 ) 


~(Py) 


where A(p y ) is the flow rule function of the inelastic phase and s are the stress deviators. The 
explicit form of the flow rule function is given by 


A(p y) = D 0 exp {-n [Z \ H l (3J^)] n } / 'W 7 


(13) 


where n = 0.5 (n + l)/n, and jf y) = yS (PY) -s (PY) is the second invariant of the average stresses in 

the inelastic subcell. Do and n are inelastic parameters, and is a state variable given for an 
isotropic hardening material by 


Z( W = zf» + (Zjp» - zf >) exp [-m wf > 1 / Zjf»] 


(14) 


where W[f Y ^ is the plastic work per unit volume. For anisotropic (directional) hardening, the 
corresponding effective state variable is given by (see Reference [2, p.220]). 


’> = ZiP« + q <fW J Z^dT + (1 - q<W) £ f Z^dt 

0 i,j=l o 


(15) 


where r,j = a,j/Voki 0 kJ . and Z(p Y ) is the time rate of change of the state variable for isotropic 
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hardening given by Equation (14). 

The meaning of the six parameters D 0 , Zo, Zj, n, m and q in Equations (14) and (15) is 
given below. 


Do = limiting strain rate in shear for large values of the second stress invariant J 2 

Zo = initial value of the hardening variable Z(t) which is related to the yield stress 

of the material in simple tension 

Zj = saturation value of the hardening variable for large values of stresses 

m = a parameter that controls the rate of work-hardening of the material 

n = a parameter that controls the rate sensitivity of the material 

q = a parameter that determines the relative weight of isotropic and directional 

hardening of the material, q = 1 provides a fully isotropic hardening and q = 0 
corresponds to fully directional hardening 

More information regarding the meaning and physical interpretation of these parameters can be 
found in Reference [2; pp. 216-224], 

The micromechanical analysis of the repeating volume element in the generalized method 
of cells proceeds in exactly the same manner as that for the original method. In this case how- 
ever, imposition of continuity of displacements and tractions between the individual subcells of 
the repeating unit cell, and between the given repeating cell and the surrounding cells, results in 
a system of equations which has a different appearance than that given by Equations (3)-(5). In 
fact, the generalization of the method of cells to an arbitrary number of subcells within the given 
repeating unit cell increases the complexity of the equations obtained by imposing the aforemen- 
tioned continuity conditions and requires the use of compact notation. The system of equations 
relating the microvariables to the uniformly applied composite strains, obtained from the imposi- 
tion of the continuity conditions, can be written in compact matrix notation as follows: 


Ae s - A P eP = Ke 


(16) 


— ~ p 

where A, A , and K are appropriate matrices of coefficients (for temperature-independent pro- 
perties, the elements of these matrices are constant); e s , e£ are the total and plastic strain vectors 
for all the subcells; and e are the average composite strains. The above equations can be solved 
in order to express the subcell strains in terms of the average strains: 
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e s = Ae + A p e£ 


(17) 


where A = A K and A p = A * A P . The matrix A, in fact, is the resulting Hill’s concentration 
matrix which relates the micro to macro variables. 

Using the definition of average stress in the repeating unit cell of the composite. 


a = 


1 Vh 1 s (P y) 
hi p 


(18) 


we can readily establish the overall constitutive equations for the multi-phase, inelastic compo- 
site in the form: 


a = B* (e - t? - a* AT) 


where 

a = the average composite stress 

e = the average composite total strain 

e* 1 = the average composite plastic strain 

B* = the effective elastic stiffness matrix of the composite 

a* = the effective CTE of the composite 

The effective stiffness matrix of the composite, B*, has the form. 


B 


* 


bn bn bj3 0 0 0 

t>i2 &22 t>23 0 0 0 

t>13 t>23 t>33 0 0 0 

0 0 0 b44 0 0 

0 0 0 0 b 55 0 

0 0 0 0 0 b 66 


(19) 


(20) 


Note that in the most general case, B* represents a composite material whose effective behavior 
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is orthotropic, in contrast to the transversely isotropic behavior given by Equation (9). 

The effective engineering constants of the composite can be determined by inverting B*. 
This inversion produces the effective compliance matrix S* whose elements are related to the 
engineering constants as follows: 

En = 1/Sn 

v 12 = — S 12/S 11 
E22 = 1/S 22 
V 2 3 = ~S 23 /S 22 

e; 3 = i/s; 3 

G23 = 1/S44 
g ; 3 = 1/S55 

G 12 = 1/S 66 


The matrix of the effective coefficients of thermal expansion is given by 

a* =diag[a*i, ot 2 2 » a 33’ (21) 

where an, a* 2 and a 33 ^ the effective CTE’s in the xj, x 2 and x 3 directions, respectively. 

— P 

The composite plastic strain vector e is determined by an incremental procedure in a step- 
wise manner. It is given in terms of the plastic strains in the individual subcells by an explicit 
expression which will not be reproduced here. Details of the GMC analysis can be found in 
Reference [3]. 
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PROGRAM DESCRIPTION 


The flow chart outlining the logical organization of the program gmc.f based on the gen- 
eralized method of cells is given on the following page. The current capabilities of the program 
are listed in Table I. 

The input data is read in from a file called gmc.data and the output is written to two files, 
namely gmc.out and gmc.plot. The input data is logically organized into four blocks, namely: 
the properties of different materials that occupy various subcells in the repeating unit cell; the 
loading applied to the composite; the internal microstructure of the repeating unit cell; and the 
geometrical details of the repeating unit cell. The output file gmc.out contains the echo of the 
input data, results for the effective thermoelastic properties, including elements of the effective 
stiffness matrix and the engineering properties, and the values of all the average stress and strain 
quantities calculated at the specified loading increments. 

The program contains three parameter variables that can be reset within the program if 
desired. The program must be recompiled if these parameters are reset by the user. The three 
parameters are: 

NBM = maximum number of subcells in the x 2 direction (set to 8) 

NGM = maximum number of subcells in the X3 direction (set to 8) 

NMM = maximum number of different materials within the repeating unit cell (set to 5 ) 


In the input file, the corresponding actual numbers for the user specified geometry are 
denoted by: 

NB = actual number of subcells in the x 2 direction 
NG = actual number of subcells in the X3 direction 
NMT = actual number of different materials 

These numbers should not exceed the corresponding maximum numbers. 

Detailed description of the input and output files, together with examples that illustrate the 
structure of these files and the capabilities of the program, are given on the following pages. 


11 


FLOW CHART FOR THE COMPUTER PROGRAM GMC.F 


READ INPUT FILE GCM.DATA 


material parameters, loading parameters, 
RVE geometry 



WRITE GCM.OUT 


effective thermoelastic coefficients 


WRITE GCM.OUT 


inelastic response at every NSTEP increment 


-*< IF NINA = 0 



WRITE GCM.PLOT 


stress-strain response at every 
NSTEP increment 


Table I. Current available capabilities within GMC. 


Type 

Description 



Constitutive models 

Elastic: isotropic, transversely isotropic materials 
Viscoplastic: Bodner-Partom unified theory for isotropic 
materials with isotropic or anisotropic hardening. 

Integration schemes 

Forward Euler integration technique 

Loadings 

Monotonic (strain control) 

9 types of uniaxial load histories (see Table II) 

Repeating unit cell geometry 

6 standard geometries (see Table III) and a user 
specified geometry 

Modelling capabilities 

monotonic thermo-mechanical response 

modelling of variable fiber shapes 

analysis of different fiber arrays 

modelling of porosities and damage 

modelling of interface regions around inclusions and fibers 

Predictive capabilities 

effective (macroscopic) elastic moduli 
effective coefficients of thermal expansion 
effective thermal conductivities 
effective stress-strain response in the inelastic region 
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Input file gmc.data 

The input file data is organized into four distinct blocks. The organization of the input data 
file and the variable names read by the program are given below. 


Block 1 - Defines constituent material parameters 


NMT - # of different materials 


Repeat 

NMT 

times 


MNAME 

EA,ET,GA 

FNA,FNT 

ALPA,ALPT 

D0,Z0,Z1 

BM,AN,Q 


- material name 

* elastic Young’s and shear moduli 

- axial and transverse Poisson’s ratios 

- axial and transverse thermal expansion coefficients 

- Bodner-Partom viscoplastic model parameters 2 

- Bodner-Partom viscoplastic model parameters 


Block 2 - Defines loading parameters 3 


5 

lines 


AMP.RATE 

LOP,NTNT,NINA 

NSTEP 

THETA 

DTEMP 


- final strain amplitude, strain rate of loading 

- load option (see Table II), number of load increments, 
actual number of load increments 

- output intervals 

- angle for LOP = 9 (skip if LOP is other than 9) 

- temperature deviation from reference temperature 


Block 3 - Defines microstructure (assigns materials to subcells in the repeating unit cell). 


NB,NG - number of subcells in the repeating unit cell in the xt and X 3 directions. 


Repeat r X...NB 

NG J 

times I 


- elements of the matrix MATNUM(NG,NB) containing material numbers 
corresponding to the material occupying the given subcell, starting with the 
subcell in the upper left comer of the repeating unit cell and progressing 
through NB columns (see Example 1). Each element of MATNUM specifies 
the material number MNAME(NM) in the given subcell (py), see Table ID 


2 Note that the material in the given subcell is perfectly elastic if DO = 0.0. 
3 See footnotes a. b. and c on following page 
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Block 4 - Defines subcell geometry 


ID - standard geometry definition (see Table III). The following are 
currently available: 1,2, 3,4,5, 6 or 10 

VF - volume fraction (needed unless ID = 10). If ID = 10, the internal cell 
geometry must be input through the specification of the dimension 
of the individual subcells. In this case, we need to specify hp and ly 
for P = 1, • • • ,Np and y = 1, • • • ,Ny, Figure 3, as follows: 



lline l\ h ... ly 


Footnotes 

a. ) 

AMP is the final value ( amplitude ) of applied (monotonic) loading. For example, if we want to run the 
program in order to generate the response of the composite in the x 2 direction, then a stress-strain curve 
d 22 “ ^22 w generated which gives the average composite stress d 22 for the given average composite 
strain e 2 2 For a final value oft 22 =0.02, say, AMP is input as 0.02. RATE is the rate of straining which is 
imposed in generating the d 22 - € 22 curve. For the rate of 0. OH sec, RATE is input as 0.01. 

b. ) 

HINT is chosen such that the incremental integration procedure for the visco-plastic model converges. 
Since the response of an inelastic composite is generated by an incremental procedure, according to 
which the strain is changed in a step-wise manner, we need to define the number of increments required to 
achieve the final value of the applied strain (i.e., to reach AMP). This number is given by NINT. In case of 
nonconvergence, the user has to increase NINT. At this stage of the program development, 
nonconvergence is ascertained by the generation of very large values for the streses and strains. An 
automatic nonconvergence criterion will be incorporated into the program in the next phase of program 
development. 


c.) 


NINA is the actual number of increments that the program will execute and may be different from NINT. If 
the effective elastic moduli and coefficients of thermal expansion are the only quantities sought by the 
user, set NINA = 0. 
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Table II. Currently available loading options. 


LOP Number 

Description 

LOP = 1 

applies £j i for a nonzero ©n while keeping all other 
stresses zero (simulates uniaxial loading in the x t 
direction) 

LOP = 2 

applies £22 for a nonzero O22 while keeping all other 
stresses zero (simulates uniaxial loading in the X2 
direction) 

LOP = 3 

applies £33 for a nonzero O33 while keeping all other 
stresses zero (simulates uniaxial loading in the X3 
direction) 

LOP = 4 

applies £23 for a nonzero 023 while keeping all other 
stresses zero (simulates transverse shear loading in the 
X2 - X3 plane) 

LOP = 5 

applies £33 for a nonzero G13 while keeping all other 
stresses zero (simulates axial shear loading in the x\ - X3 
plane) 

LOP = 6 

applies £12 for a nonzero On while keeping all other 
stresses zero (simulates axial shear loading in the xj - X2 
plane) 

LOP = 7 

applies £22 = -£33 for nonzero G22 and 033 while keeping 
all other stresses zero (simulates pure shear in the X2 - X3 
plane at 45 0 to the applied loading) 

LOP = 8 

applies AT (temperature deviation from a reference 
temperature) while keeping all composite stresses Ojj zero. 

LOP = 9 

applies £ xx to a unidirectional lamina in which the fibers 
are inclined at an angle 9 with respect to the loading 
direction for a nonzero o xx , while keeping all other stresses 
zero (simulates axial loading of an off-axis lamina). 
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Table III. Standard repeating unit cell geometry definitions (see Appendix I for a more detailed 
explanation). 


ID Number 

Geometry configuration 

Material assignment matrix 4 

ID = 1 

"square" fiber in a square array 

2 2 


(see Figure 6) 

1 2 

ID = 2 

"square" fiber in a triangular array 

2 2 12 


(see Figure 7) 

2 2 2 2 


VF < 0.288675 (see Figure 8a) 

12 2 2 
2 2 2 2 



12 12 


VF > 0.288675 (see Figure 8b) 

12 2 2 
12 12 
2 2 12 

ID = 3 

"square" fiber in a square diagonal array 

2 2 12 


(see Figures 9 and 10) 

2 2 2 2 
12 2 2 
2 2 2 2 

ID = 4 

cross-shaped fiber in a square array 

2 2 2 2 


(see Figures 1 1 and 12) 

2 12 2 
1112 
2 12 2 

ID = 5 

coated, "square" fiber in a square array 

2 2 2 2 


(see Figures 13 and 14) 

3 3 3 2 
3 13 2 
3 3 3 2 

ID = 6 

circular fiber in a square array 

2 2 2 2 2 2 2 


(see Figure 15) 

2 2 2 1 2 2 2 
2 2 1112 2 
2 111112 
2 2 1112 2 
2 2 2 1 2 2 2 
2 2 2 2 2 2 2 

ID= 10 

user specified fiber shape and geometry 

any assignment 


4 Material 1 - fiber, material 2 - matrix, material 3 - coating 
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Example 1 

This example illustrates the construction of the matrix MATNUM(NB,NG) that stores 
information on the types of materials that occupy various subcells. Each of the elements of 
MATNUM specifies the material number MNAME(NM) in the given subcell (|3y). Let NMT = 
2, NB = 3, NG = 4. For the specified NB and NG, the repeating unit cell has the form shown in 
Figure 4. In order to show the form of MATNUM(NB,NG), let us assign in this example 
’Material 1’ to the subcells (1,1), (2,2) and (1,4) while allowing the remaining subcells to be 
occupied by ’Material 2’. The resulting assignment is depicted in Figure 5. 

Example 2 

This example shows how to construct an input file in order to generate the transverse 
stress-strain response curve, 022 versus £ 22 0- e - LOP = 2), of a boron/aluminum composite using 
a repeating unit cell that consists of 4 subcells (i.e. NB = 2 and NG = 2, ID = 1). The properties 
of the elastic boron fibers, ’Material 1’, and the elastic-perfectly plastic aluminum matrix, 
’Material 2’, are given as talcen from Reference [2], p. 248. They are reproduced in the table 
below. 


Table IV. Material parameters of boron fibers and aluminum matrix 


Material 

E 

(GPa) 

V 

Do 1 

(sec) 

Zo 

(MPa) 

Zi 

(MPA) 

m n 

Boron 

Aluminum 5 

413.70 

55.16 

0.20 

0.30 

10 -4 

103.42 

103.42 

10 


The input file is given in Appendix n. Note that AMP = 0.005 and RATE = 0.01/sec. The 
number of increments required to generate the response curve is NINT = 2000. This number 
ensures that the resulting integration time increment is sufficiently small to attain convergence. 
The output is written to the output file gmc.out every 20 increments (since NSTEP = 20). In 
addition, DTEMP = 0 (i.e. AT = 0.0) and VF=0.5. 


5 Note that since Z 0 = Z t . the material is elastic, prefectly plastic and thus m is irrelevant (i.e., can be anything). 
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Example 3 

For the same boron/aluminum system of Example 2 we generate the transverse stress-strain 
response curve, G 2 2 versus £ 22 , using 4x4 subcells in the repeating unit cell. As the fundamen- 
tal configuration and geometry of the repeating unit cell shown in Figure 16 produces the same 
doubly periodic array as in Example 2, the resulting stress-strain response curve should be ident- 
ical. Here, NB = 4, NG = 4 and VF = 0.5. 

Since the geometry in this case is supplied by the user, ID = 10 and so the dimensions hp 
and 1 Y of the subcells ((3y) must be specified. In order to generate VF = 0.5, we take 


hi + lb — 0.5 ; li +I 2 = 0-5 

hi = h 3 ; li = I 3 


h 2 = h 4 ; 1 2 = I 4 

h] = lj : h 2 = 1 2 


Since VF = 


4xhi 


-> hi = — = 0.3535. The subcell dimensions are input in the following 

O 


manner 


h 4 

h 3 

h 2 

hi 


or 

0.1464 

0.3535 

0.1464 

0.3535 

0.3535 0.1464 0.3535 0.1464 

The complete form of the input file for the present example is given in Appendix III. 


19 


Output file gmc.out 

The first section of this file contains an "echo" of the input data that was read by the pro- 
gram from the input file gmc.data. The next section of the output file contains information about 
the subcell identification, material assignment in each subcell, volume of each subcell and the 
total volume of the repeating unit cell. It is a concise summary about the internal arrangement of 
the individual subcells. 

The next section presents the effective stiffness matrix B* of the composite. This is fol- 
lowed by the effective engineering elastic moduli: 

E11S ~> En 
N12S ~> v *2 
E22S ~> Eo 2 
N23S ~> V 23 
E33S ~> E 33 
G23S ~> G 23 
G13S ~> G 13 
G12S ~> Gl 2 


The effective coefficients of thermal expansion are printed as follows: 

ALPHA1 IS -> a„ 

ALPHA22S -> a 22 
ALPHA33S -> 033 


Note 1 : If the effective elastic moduli and coefficients of thermal expansion are the only quanti- 
ties sought by the user, the parameter NINA = 0 in the input file as mentioned previously. In this 
case, the program will not generate the inelastic response and the output file will terminate at 
this stage. 

Note 2 : If the effective thermal conductivities of the composite k*i , k 22 and k 33 , are sought, the 
user should set NINA = 0, and the following replacement should be made in the input file: 


EA(NM) — > k A for each material type 
ET(NM) — > kj for each material type 
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FNA(NM) = FNT(NM) = 0 
GA(NM) is arbitrary 

where k A and kx are the axial and transverse thermal conductivities of a transversely isotropic 
matarial which occupies a given subcell. In the output file, k^, V-n and ^33 316 identified by 
Ejj, E 22 and E 33 , respectively. 

The next section presents the inelastic response when the parameter NINA is greater than 
zero. The program prints the following results at the first increment ( INC = 7) and subsequently 
at every NSTEP increments according to the format: 

INC = increment number 


IJ component Total strains Stresses Plastic strains Thermal strains 

11 

22 

33 

23 

13 

12 


The program is terminated when the increment number INC reaches NINA. A sample 
gmc.out output file is presented in Appendix IV generated with the input file gmc.data of 
Example 3. 


21 


Output file gmc.plot 

For the parameter NINA greater or equal to 1 , the program writes to this file the specified 
stress-strain response dictated by the loading option LOP , starting at the first increment INC = 1, 
and at each subsequent NSTEP increment. The loading option specifies the stress-strain response 


written to this 

file according 

to; 

LOP = 1 

£ 11 

611 


LOP = 2 

£22 

o 2 2 


LOP = 3 

£33 

633 


LOP = 4 

^23 

(J23 


LOP = 5 

£13 

^13 


LOP = 6 

£12 

Oi 2 

y (^22 

LOP = 7 

—(C22 - £33) 

LOP = 8 

AT 

£jl £22 

£33 

LOP = 9 

£xx 

°xx 
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APPENDIX I 


A more detailed explanation of the standard repeating unit cell geometry definitions sum- 
marized in Table II are given in this appendix. 


ID = 1 : corresponds to the original method of cells (NB = NG = 2) in which the fiber 
occupies subcell (1,1), with the other three subcells occupied by the matrix as shown in 
Figure 6. Let VF denote the fiber volume fraction. The dimension of the repeating unit 
cell are scaled to lxl so that hj + h2 = li + 12 = 1. Thus, for a "square" fiber, 
speci ficati on of VF is sufficient to define the dimension hp and ly, i.e., 
hj = VVF , I12 = 1 - hi , lj = h] and U = h2- The material assignment of the fiber 
(material 1 ) and matrix (material 2) is as follows: 


2 2 
1 2 


ID = 2 : gives a triangular packing of "square" fibers shown in Figure 7. In this 
arrangement, NB = NG = 4, and for a given fiber volume fraction VF the dimensions hp and 
ly are automatically determined internally. For this packing arrangement, it is necessary to 
differentiate between two cases which is done internally. For VF < 0.288675, the geometry 
of the repeating unit cell is given in Figure 8a, where 

hi = b-| 

ho = a 
h3=h, 
h 4 = a 

and 

h =a 

1 2 =c 

1 3 =a 

14 =c 


where 


a= Vi^ VF 

c= VF- a 

The material assignment of the fiber (material 1) and matrix (Material 2) in this case is: 
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2 2 12 
2 2 2 2 
12 2 2 
2 2 2 2 


For VF > 0. 288675, the geometry of the repeating unit cell is given in Figure 8b, where 

hj =2b 

h 2 = f -b 

h 3 =h, 
h 4 = h 3 

and 

h=a 


The material assignment of the fiber (material 1 ) and matrix ( material 2) in this situation is: 


ID — 3 : specifies a square diagonal packing illustrated in Figure 9. In this case NB = NG = 
4. The subcell dimensions hp and ly are determined internally in the program for the 
specified fiber volume fraction VF. For this packing arrangement, VF cannot exceed 0.5. 
The dimension of the repeating unit cell for this packing arrangement shown in Figure 10, 
where 



12 12 
12 2 2 
12 12 
2 2 12 



h 3 =h 1 
h 4 — h 4 


and 
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l,=h 2 

h =hj 

b =h 

1 4 = 1 2 

The material assignment of the fiber (material 1) and matrix (material 2) for this case is: 


2 2 12 
2 2 2 2 
12 2 2 
2 2 2 2 


ID - 4 : models a cross-shaped fiber shown in Figure 11. It requires that NB = NG = 4, and 
the parameter "a” shown in the figure must be provided by the user. Thus in this case the 
parameter "a" should be read by the statement READ A. For this arrangement, the 
geometric details o f the repeating unit cell are shown in Figure 12. Defining 
x = -2a + v4a 2 + VF , the dimensions hp and ly are given by 

hj = a 
h 2 = x 
h 3 = a 

h 4 = 1 - 2a - x 
and 

li=hi 

1 2 = h 2 

13 =h 3 

1 4 = h 4 

The material assignment of the fiber (material 1 ) and matrix (material 2) for this case is: 


2 

2 

2 

2 

2 

1 

2 

2 

1 

1 

1 

2 

2 

1 

2 

2 


ID - 5 : models a coated, "square" fiber shown in Figure 13. Let VF and VC denote, 
respectively, the fiber volume fraction and the coating volume fraction. This situation 
requires that NB = NG = 4. Note that in this case the parameter VC is read by the statement 
READ VC. The dimensions of the repeating unit cell for this arrangement are given in 
Figure 14, where 
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hj = — [-VVF + VVF + VC] 
h 2 = VVF 

h3 = hi 

h4 = 1 - [hi + h 2 + 113] 


and 

11 =hi 

1 2 = h 2 
h - h 3 

I4 = h4 


The material assignment of the fiber (material 1), matrix (material 2), and coating (material 
3) for this case is: 


2 2 2 2 

3 3 3 2 

3 13 2 

3 3 3 2 


ID = 6 : this case approximates a circular fiber using Np = 7 and Ny = 7 subcells as shown 
in Figure 15. By imposing the condition that Vf = TtR'/h", where R is the fiber radius, it 
follows that 

h 2 = h 3 = h 5 = h 6 = R 
h 4 = 4h 2 

h, = h 7 = [h - (2h 2 + 2h 3 + h 4 )]/2 and 1, = hjor i = 1, • • • ,7 

The material assignment of the fiber (material 1 ) and matrix (material 2) for this case is: 


2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

1 

2 

2 

2 

2 

2 

1 

1 

1 

2 

2 

2 

1 

1 

1 

1 

1 

2 

2 

2 

1 

1 

1 

2 

2 

2 

2 

2 

1 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 


ID = 10 : enables the user to specify his own internal configuration by specifying the 
dimensions hp and ly. 
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APPENDIX II 


The input file for the case described in Example 2 is given below. 


2 

'Material 1' 




413 .70E+09 

413.70E+09 

172 . 30E+09 



0 .20E+00 

0 .20E+00 




1 . 00E-06 

1 . 00E-06 




0 . 00E+00 

1 . 00E+00 

1 .00E+00 



1 . OOE+OO 

1 .00E+00 

1 . 00E+00 

Block 

1 - 

'Material 2 

t 




55 . 16E+09 

55.16E+09 

21 . 20E+09 



0.30E+00 

0 . 30E+00 




1 . 00E-06 

1.00E-06 




» 1.00E+04 

103. 42E+0 6 

103 . 42E+06 



1 .00E+00 

10 . 00E+00 

1 . 00E+00 



f 0.005 0. 

01 




J 2 2000 

2000 


Block 

2 

1 20 





(. 0.0 





( 2 2 



Block 

3 

1 2 2 





l 1 2 





f 1 



Block 

4 

1 0.5 
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APPENDIX III 


The input file 


Block 1 


Block 2 


Block 3 


Block 4 


for the case described in Example 3 is given below. 


'Material 1' 


413 . 70E+09 

413.7 0E+0 9 

172.30E+09 

0 . 20E+00 

0 .20E+00 


1.00E-06 

1.00E-06 


0 . 00E+00 

1 . 00E+00 

1 . OOFi-OO 

1.00E+00 

1 .00E+00 

1.00E+00 

'Material 2 

r 


55.16E+09 

55 . 1 60E+0 9 

21 . 20E+09 

0.30E+00 

0 . 30E+00 


1 . 00E-06 

1 . 00E-06 


1 . 00E+04 

103 . 42E+0 6 

103 . 42E+06 

1 . 00E+00 

10 . 00E+00 

1 . 00E+00 

0.005 0. 

01 


2 2000 

2000 



20 

0.0 


4 
2 
1 
2 
1 
10 

0 .1464 
,3535 
, 1464 
,3535 
,3535 


0.1464 0.3535 0.1464 


APPENDIX IV 


The output file gmc.out generated by the input file presented in Example 3 is given below. 


**************************************************** 

** GENERALIZED CELLS MODEL ** 

★ * * * 


★ * 
* * 
★ * 
* * 
★ * 
★ * 
★ * 


G C M ** 

★ * 

Determination of effective elastic, thermal ** 

and inelastic response of unidirectional, ** 

multi-phase composites using generalized ** 

version of the method of cells ** 

★ * 


** Programmed by ** 
** Jacob Aboudi & Moshe Paley ** 
** NASA Lewis version - July 1992 ** 
**************************************************** 


***** input data echo ***** 


MATERIAL SPECIFICATION 


Number of Materials (NMT ) = 2 


Material 1 


Thermo-elastic constants 


EA = 0.414E+12 

ET 

= 

0 . 414E + 12 

GA = 

0 . 172E+12 

NUA = 0 . 200E+00 

NUT 

= 

0 .200E+00 



ALFA = 0 . 100E-05 

ALFT 

= 

0 . 100E-05 



Inelastic parameters 





DO = 0 . 000E+00 

ZO 

= 

0 . 100E+01 

Z1 = 

0 . 100E+01 

M = 0 . 100E+01 

N 

= 

0 . 100E+01 

Q = 

0 . 100E+01 

Material 2 






Thermo-elastic constants 





EA = 0.552E+11 

ET 

= 

0 . 552E+11 

GA = 

0 .212E+11 

NUA = 0 . 300E+00 

NUT 

= 

0 . 300E+00 



ALFA = 0 . 100E-05 

ALFT 

= 

0 . 100E-05 



Inelastic parameters 





DO = 0 . 100E+05 

ZO 

= 

0 . 103E+09 

Z1 = 

0 . 103E+09 

M = 0.100E+01 

N 

= 

0 . 100E+02 

Q = 

0 . 100E+01 
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LOADING SPECIFICATION 


Amplitude = 0.500E-02 Rate = 0.100E-01 

Loading option = 2 

Number of integration increments = 2000 

Actual number of increments = 2000 

Number of print steps = 20 

Temperature change = 0.000E+00 


RVE SPECIFICATION 


Number of subcells in the x-2 direction = 4 

Number of subcells in the x-3 direction = 4 


Subcell material assignment 


2 2 2 2 
12 12 
2 2 2 2 
12 12 


Subcell dimension assignment 
ID = 10 


H 4 = 0.14 6E + 0 0 
H 3 = 0 . 354E+00 
H 2 = 0.14 6E + 0 0 
HI = 0 . 354E+00 

L(1,...,NG) = 0 . 354E+00 0.146E+00 0.354E+00 0.146E+00 


***** OUTPUT RESULTS ***** 


SUBCELL IDENTIFICATION, VOLUME AND MATERIAL ARRANGEMENT 


(BETA 

f 

GAMMA) 

SUBCELL # 

SUBCELL MATERIAL 

SUBCELL VOLUME 

1 

r 

1 

1 

1 

0 . 125E+00 

1 

r 

2 

2 

2 

0 .518E-01 

1 

r 

3 

3 

1 

0 . 125E+00 

1 

r 

4 

4 

2 

0 . 518E-01 

2 

r 

1 

5 

2 

0.518E-01 

2 

r 

2 

6 

2 

0 .214E-C1 

2 

r 

3 

7 

2 

0 . 518E-01 

2 

f 

4 

8 

2 

0 .214E-01 

3 

t 

1 

9 

1 

0 . 125E+00 

3 

t 

2 

10 

2 

0.518E-01 

3 

! 

3 

11 

1 

0 . 125E+00 

3 

t 

4 

12 

2 

0.518E-01 

4 

f 

1 

13 

2 

0 . 518E-01 

4 

t 

2 

14 

2 

0 .214E-01 
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2 0 . 518E-01 

2 0 . 214E-01 


TOTAL VOLUME = 0.100E+01 


Effective stiffness matrix 


0 . 258E+12 0 . 473E+11 0.473E+11 

0 . 473E+11 0.14 9E+12 0.446E+11 

0 . 473E+11 0.446E+11 0.149E+12 

0 . 378E+11 

0 . 457E+11 

0 . 457E+11 


Effective engineering moduli 


El 1S= 0 . 235E+12 
N12S= 0 . 245E+00 
E22S= 0 . 131E+12 
N23S= 0 . 257E+00 
E33S= 0 . 131E+12 
G23S= 0 . 378E+11 
G13S= 0 . 457E+11 
G12S= 0 . 457E+1 1 


4.3 15 

4.4 16 


Effective thermal expansion coefficients 


ALPHA1 1S= 
ALPHA22S= 
ALPHA33S= 


0 . 100E-05 
0 . 100E-05 
0 . 100E-05 


EFFECTIVE STRESS-STRAIN RESPONSE 


INC = 1 


IJ 

TOTAL STRAINS 

STRESSES 

PLASTIC STRAINS 

THERMAL STRAINS 

11 

-0 . 341E-06 

0 . 195E-02 

0 . 000E+00 

0 . 000E+00 

22 

0 .250E-05 

0 . 327E+0 6 

0 . 000E+00 

0 . 000E+00 

33 

-0.642E-06 

0 . 391E-01 

0.000E+00 

0 . 000E+00 

23 

0 . 000E+00 

0 . 000E+00 

0.000E+00 

0.000E+00 

13 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

12 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 


INC 

= 20 




IJ 

TOTAL STRAINS 

STRESSES 

PLASTIC STRAINS 

THERMAL STR 
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11 

-0 . 682E-05 

0 . 625E-01 

-0 . 504E-21 

0. 000E+00 

22 

0.500E-04 

0 . 654E+07 

0 . 328E-20 

0 . 000E+00 

33 

-0 . 128E-04 

0 . 500E+00 

-0 .232E-20 

0. 000E+00 

23 

0 . 000E+00 

0 .000E+00 

0 . 000E+00 

0. 000E+00 

13 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

0. 000E+00 

12 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

0 .000E+00 

INC 

= 40 




IJ 

TOTAL STRAINS 

STRESSES 

PLASTIC STRAINS 

THERMAL STR 

11 

-0 . 136E-04 

0 . 125E+00 

-0 . 103E-20 

0 .000E+00 

22 

0 . 100E-03 

0 . 131E+08 

0 . 673E-20 

0 . 000E+00 

33 

-0 .257E-04 

0 . 100E+01 

-0 . 476E-20 

0 . 000E+00 

23 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

0 . 000E+00 

13 

0 . 000E+00 

0 . 000E+00 

0 . OOOE+OO 

0 . 000E+00 

12 

0 . 000E+00 

0 . 000E+00 

0 .000E+00 

0 . 000E+00 


INC 

= 1980 




IJ 

TOTAL STRAINS 

STRESSES 

PLASTIC STRAINS 

THERMAL STRAINS 

11 

-0.316E-03 

-0 . 400E+01 

-0 . 121E-03 

0 . OOOE+OO 

22 

0 . 4 95E-02 

0 . 187E+09 

0.352E-02 

0 .OOOE+OO 

33 

-0.355E-02 

0 . 240E+02 

-0 . 318E-02 

0 . OOOE+OO 

23 

0 . 000E+00 

0. 000E+00 

0 . OOOE+OO 

0 . OOOE+OO 

13 

0 .000E+00 

0 .000E+00 

0 . OOOE+OO 

0 . OOOE+OO 

12 

0. 000E+00 

0. 000E+00 

0. OOOE+OO 

0 . OOOE+OO 

INC 

= 2000 




IJ 

TOTAL STRAINS 

STRESSES 

PLASTIC STRAINS 

THERMAL STRAINS 

11 

-0 . 31 6E-03 

-0 . 600E+01 

-0 . 121E-03 

0. OOOE+OO 

22 

0 . 500E-02 

0 . 187E+09 

0 . 3S7E-02 

0 . OOOE+OO 

33 

-0 . 3 60E-02 

0 .000E+00 

-0.323E-02 

0 . OOOE+OO 

23 

0. 000E+00 

0 .000E+00 

0 . OOOE+OO 

0 . OOOE+OO 

13 

0 .000E + 00 

0 .OOOE+OO 

0 . OOOE+OO 

0 . OOOE+OO 

12 

0 .000E+00 

0 . OOOE+OO 

0 . OOOE+OO 

0 .OOOE+OO 
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a) Doubly periodic array 


b) Representative unit cell 


Figure 1 . Doubly-periodic array and a repeating unit cell for the original method of cells. 
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Figure 2. Unit periodic geometry for the generalized method of cells (GMC). 
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Figure 3. A repeating volume element of GMC. 
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Figure 4. A repeating volume element with NB = 3 and NG = 4 (Example 1). 
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Figure 5. Material assignment for Figure 4. 
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Figure 6. Repeating volume element for the case ED = 1, with the fiber assigned 
to the subcell (1,1) and the matrix to the remaining subcells. 



Figure 7. Fiber arrangement for triangular packing (ID = 2). 
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h l 2 I 3 l 4 


a.) 



*1 l 2 h Z 4 


b.) 

Figure 8. Fiber and matrix locations in a repeating unit cell for triangular 
packing (ID = 2) when a.) v f < 0.288675 and b.) Vf > 0.288675. 


Figure 9. Fiber arrangement for square diagonal packing (ID = 3). 



■< >■ 

h l 2 l 3 U 


Figure 10. Fiber and matrix locations in a repeating unit cell for square diagonal 

packing (ID = 3). 



Figure 11. A cross-shaped fiber (ID = 4). 



/l / 2 /3 I 4 


Figure 12. Fiber and matrix locations in a repeating unit cell for a cross-shaped 

fiber (ID = 4). 
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Figure 13. A fiber with a coating (ID = 5). 



Coating 

Fiber 

Matrix 
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Figure 14. Fiber, matrix and coating locations in a repeating unit cell (ID = 5). 
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Figure 15. Repeating volume element for a circular fiber approximation (ID = 6). 



Figure 16. Material assignment for Example 3. 
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